library(Matrix)	
library(eliot)
setwd("~/Dropbox/Shared Folders/Ben Lauderdale - LDA+IRT/CodeForPaper/")
#setwd("~/Dropbox/Research/Project - LDA+IRT (TC)/CodeForPaper/")	
#source("LDA+IRTFunctions.R")
load("LDA+IRTPaperData.RData")


## Plot DIC as a Function of Topic Count ##

SavedTopicCountRuns <- list.files("~/Dropbox/Research/Project - LDA+IRT (Saved MCMC)/")
SavedTopicCountRuns <- SavedTopicCountRuns[grep("2000-4000",SavedTopicCountRuns)]
Nruns <- length(SavedTopicCountRuns)
DICmat <- matrix(NA,Nruns,2)
for (i in 1:Nruns){
	load(paste("~/Dropbox/Research/Project - LDA+IRT (Saved MCMC)/",SavedTopicCountRuns[i],sep=""))
	temptopics <- substr(SavedTopicCountRuns[i],13,14)
	if (substr(temptopics,2,2) == "-") temptopics <- substr(temptopics,1,1)
	DICmat[i,1] <- temptopics
	DICmat[i,2] <- LDAIRT.out$dic.irt	
	}

pdf(file="Plots/OptimalTopicCount2.pdf",6,4)
par(mfrow=c(1,1),mar=c(5.1,4.1,4.1,2.1))
plot(DICmat[,1],DICmat[,2],xlab="Topic Count",ylab="DIC",main="DIC for Roll Call Votes by Topic Count",ylim=c(27500,30750),axes=FALSE,type="h")
axis(1,at=seq(0,40,5))
axis(2)
box()
dev.off()


## Load Best MCMC Chain and Create Summaries ## 

#load(file="~/Dropbox/Research/Project - LDA+IRT (Saved MCMC)/SavedLDAIRT-24-5000-10000.RData")
load(file="OptimalTopicCountSearch/SavedLDAIRT-24-5000-5000.RData")
LDAIRT.summary <- summary.LDAIRT(LDAIRT.out)
attach(LDAIRT.summary)
theta.se <- apply(LDAIRT.out$ideal.chain,c(1,2),sd)

SpaethCategoryNames <- c("Criminal Procedure","Civil Rights","First Amendment","Due Process","Privacy","Attorneys","Unions","Economic Activity","Judicial Power","Federalism","Interstate Relations","Federal Taxation","Miscellaneous")
LDATopicLabels <- rep("",Ntopics)
for (k in 1:Ntopics) LDATopicLabels[k] <- paste((WordSparseMatrix$dimnames$Terms[TopWordsByTopic[,k]])[1:3],collapse=", ")
LDATopicTable <- matrix("",10,Ntopics)
for (k in 1:Ntopics) LDATopicTable[,k] <- WordSparseMatrix$dimnames$Terms[TopWordsByTopic[1:10,k]]

## Save Distinctive Word Table ##

write.table(t(LDATopicTable),file="Plots/TopWordTable.tex",sep="&",quote=FALSE,row.names=FALSE,col.names=FALSE)

## Generate Example Plot ##

BGCourt <- which(is.na(VoteMatrix[4242,]) == FALSE)
BGNames <- colnames((lambda.est[,4242] %*% theta.est))[BGCourt]
BGtheta <- (lambda.est[,4242] %*% theta.est)[BGCourt]
names(BGtheta) <- BGNames
BGlambda <- lambda.est[,4242]









## Plot Heat Map of Topics By Spaeth Code and Spaeth Code by LDA Topics ##

pdf(file="Plots/TopicHeatMaps1.pdf",8,8)
par(mfrow=c(1,1),mar=c(12,10,4,4))
SpaethTopicWeightMatrix <- matrix(NA,Ntopics,13)
for (k in 1:Ntopics){
	SpaethTopicWeightMatrix[k,] <- lm(lambda.est[k,]~as.factor(VoteDetails$IssueArea)-1)$coef
}
# plot(0,0,type="l",xlim=c(0,13),ylim=c(0,Ntopics),xlab="",ylab="",axes=FALSE,main="LDA Topic Weight by Spaeth Issue Area")
# axis(1,at=seq(0.5,12.5,1),labels=SpaethCategoryNames,las=2)
# axis(2,at=seq(0.5,Ntopics-0.5,1),labels=LDATopicLabels,las=2,cex.axis=0.65)
# for (k in 1:Ntopics){
# 	for (sp in 1:13){
# 		rect(sp-1,k-1,sp,k,col=grey(1-SpaethTopicWeightMatrix[k,sp]))
# 	}	
# }

TopicSpaethWeightMatrix <- matrix(0,13,Ntopics)
for (j in 1:Nvotes){
	for (k in 1:Ntopics){
		TopicSpaethWeightMatrix[VoteDetails$IssueArea[j],k] <- TopicSpaethWeightMatrix[VoteDetails$IssueArea[j],k] + lambda.est[k,j]
	}
}

TopicSpaethWeightMatrix <- TopicSpaethWeightMatrix/t(matrix(colSums(TopicSpaethWeightMatrix),Ntopics,13))

plot(0,0,type="l",ylim=c(0,13),xlim=c(0,Ntopics),xlab="",ylab="",axes=FALSE,main="Spaeth Issue Area by LDA Topic")
axis(2,at=seq(0.5,12.5,1),labels=SpaethCategoryNames,las=2)
# axis(1,at=seq(0.5,Ntopics-0.5,1),labels=LDATopicLabels,las=2,cex.axis=0.65)
axis(1,at=seq(0.5,Ntopics-0.5,1),labels=FALSE,las=2,cex.axis=0.65)
for (k in 1:Ntopics){
	 for (sp in 1:13){
		 rect(k-1,sp-1,k,sp,col=grey(1-TopicSpaethWeightMatrix[sp,k]))
	 }	
 }
text(x = seq(1,Ntopics,1), par("usr")[3] - 0.5, labels = LDATopicLabels, srt = 45, pos = 2, xpd = TRUE)

dev.off()	


pdf(file="Plots/SpaethTopicHeatMap.pdf",8,8)
par(mfrow=c(1,1),mar=c(10,10,4,4))
SpaethTopicWeightMatrix <- matrix(NA,Ntopics,13)
for (k in 1:Ntopics){
	SpaethTopicWeightMatrix[k,] <- lm(lambda.est[k,]~as.factor(VoteDetails$IssueArea)-1)$coef
}
plot(0,0,type="l",xlim=c(0,13),ylim=c(0,Ntopics),xlab="",ylab="",axes=FALSE,main="LDA Topic Weight by Spaeth Issue Area")
axis(1,at=seq(0.5,12.5,1),labels=SpaethCategoryNames,las=2)
axis(2,at=seq(0.5,Ntopics-0.5,1),labels=LDATopicLabels,las=2,cex.axis=0.65)
for (k in 1:Ntopics){
	for (sp in 1:13){
		rect(sp-1,k-1,sp,k,col=grey(1-SpaethTopicWeightMatrix[k,sp]))
	}	
}

dev.off()

pdf(file="Plots/TopicSpaethHeatMap.pdf",12,8)
par(mfrow=c(1,1),mar=c(10,10,4,4))

TopicSpaethWeightMatrix <- matrix(0,13,Ntopics)
for (j in 1:Nvotes){
	for (k in 1:Ntopics){
		TopicSpaethWeightMatrix[VoteDetails$IssueArea[j],k] <- TopicSpaethWeightMatrix[VoteDetails$IssueArea[j],k] + lambda.est[k,j]
	}
}

TopicSpaethWeightMatrix <- TopicSpaethWeightMatrix/t(matrix(colSums(TopicSpaethWeightMatrix),Ntopics,13))

plot(0,0,type="l",ylim=c(0,13),xlim=c(0,Ntopics),xlab="",ylab="",axes=FALSE,main="Spaeth Issue Area by LDA Topic")
axis(2,at=seq(0.5,12.5,1),labels=SpaethCategoryNames,las=2,cex=1.5)
axis(1,at=seq(0.5,Ntopics-0.5,1),labels=FALSE)
text(x = seq(1,Ntopics,1), par("usr")[3] - 0.5, labels = LDATopicLabels, srt = 45, pos = 2, xpd = TRUE)
for (k in 1:Ntopics){
	 for (sp in 1:13){
		 rect(k-1,sp-1,k,sp,col=grey(1-TopicSpaethWeightMatrix[sp,k]))
	 }	
 }
dev.off()	







library(gplots)
library(ggplot2)
library(lattice)
library(RColorBrewer)
low <- col2rgb("white")/255
high <- col2rgb("black")/255
r <- seq(low[1], high[1], len=50)
g <- seq(low[2], high[2], len=50)
b <- seq(low[3], high[3], len=50)
pallette <- rgb(r, g, b)

pdf('Plots/SpaethTopicHeatmap2.pdf', 6, 6)
heatmap.2(SpaethTopicWeightMatrix
	, col=function (n) colorpanel(n, "white", "grey", "black")
	, distfun = function(SpaethTopicWeightMatrix)dist(SpaethTopicWeightMatrix)
	, hclustfun = function(SpaethTopicWeightMatrix) hclust(SpaethTopicWeightMatrix, method='average')
	, labCol=SpaethCategoryNames
	, labRow=LDATopicLabels
	, dendrogram='none'
	, cexCol = 0.5, cexRow = 0.75, scale = 'column', margin=c(4,10)
	, trace='none'
	, key=FALSE
	, keysize=.3
	, margins = c(5,10)
	, sepwidth = c(0,0)
	)
dev.off()

pdf('Plots/TopicSpaethHeatmap2.pdf', 6, 6)
heatmap.2(TopicSpaethWeightMatrix
	, col=function (n) colorpanel(n, "white", "grey", "black")
	, distfun = function(TopicSpaethWeightMatrix)dist(TopicSpaethWeightMatrix)
	, hclustfun = function(TopicSpaethWeightMatrix) hclust(TopicSpaethWeightMatrix, method='average')
	, labRow=SpaethCategoryNames
	, labCol=LDATopicLabels
	, dendrogram='none'
	, cexCol = 0.5, cexRow = 0.75, scale = 'column', margin=c(4,10)
	, trace='none'
	, key=FALSE
	, keysize=.3
	, margins = c(7,7)
	, sepwidth = c(0,0)
	)
dev.off()

## Plot of Issue Frequency Over Time ##


IssueSubsetForPresentation <- c(10,16,21)


pdf(file="Plots/TimeTrendsInTopics.pdf",18,12)
par(mfrow=c(4,6),mar=c(5.1,4.1,4.1,2.1))
for (k in 1:Ntopics){
	loess.out <- loess(LDAIRT.out$topic.chain[k,,1]~VoteDetails$Term,span=0.5)
	plot(loess.out$x,loess.out$fitted,type="l",ylim=c(0,0.2),col="black",lwd=2,main=LDATopicLabels[k],xlab="Year",ylab="Mean Case Fraction")
	}
dev.off()	

pdf(file="Plots/TimeTrendsInTopicsS.pdf",9,3)
par(mfrow=c(1,3),mar=c(5.1,4.1,4.1,2.1))
for (k in IssueSubsetForPresentation){
	loess.out <- loess(LDAIRT.out$topic.chain[k,,1]~VoteDetails$Term,span=0.5)
	plot(loess.out$x,loess.out$fitted,type="l",ylim=c(0,0.2),col="black",lwd=2,main=LDATopicLabels[k],xlab="Year",ylab="Mean Case Fraction")
	}
dev.off()	

pdf(file="Plots/TimeTrendsInTopics2.pdf",18,12)
par(mfrow=c(4,6), mar=c(1,1,3,1), oma=c(4,4,1,1))
for (k in 1:Ntopics){
	loess.out <- loess(LDAIRT.out$topic.chain[k,,1]~VoteDetails$Term,span=0.5)
	plot(loess.out$x,loess.out$fitted,type="n",ylim=c(0,0.2),col="black",lwd=2,main='',xlab="Year",ylab="Mean Case Fraction",xaxs='i',yaxs='i')
	title(main=list(LDATopicLabels[k], cex=1.6))
	polygon(c(1948,1953,1953,1948),c(0,0,.2,.2), col='grey',border=NA)
	polygon(c(1969,1986,1986,1969),c(0,0,.2,.2), col='grey',border=NA)
	lines(loess.out$x,loess.out$fitted, lwd=2)
	}
mtext(side=1,'Year',outer=TRUE, line=2)
mtext(side=2,'Mean Case Faction',outer=TRUE, line=2)
dev.off()

## trying to make a boxplot
library(beanplot)
pdf(file="Plots/TimeTrendsInTopics3.pdf",18,12)
par(mfrow=c(4,6), mar=c(1,1,3,1), oma=c(4,4,1,1))
for (k in 1:Ntopics){
  boxplot(LDAIRT.out$topic.chain[k,,1] ~ VoteDetails$Term
           ,main='',xlab="Year",ylab="Case Fraction"#,xaxs='i',yaxs='i'
           ,do.out=FALSE, do.conf=FALSE, outline=TRUE)
  title(main=list(LDATopicLabels[k], cex=1.6))
#   polygon(c(1948,1953,1953,1948),c(0,0,.2,.2), col='grey',border=NA)
#   polygon(c(1969,1986,1986,1969),c(0,0,.2,.2), col='grey',border=NA)
}
mtext(side=1,'Year',outer=TRUE, line=2)
mtext(side=2,'Case Faction',outer=TRUE, line=2)
dev.off()

VoteDecade <- ifelse(VoteDetails$Term < 1960, 1950,
                     ifelse(VoteDetails$Term>1959 & VoteDetails$Term<1970, 1960,
                            ifelse(VoteDetails$Term>1969 & VoteDetails$Term<1980, 1970,
                                   ifelse(VoteDetails$Term>1979 & VoteDetails$Term<1990, 1980,
                                          ifelse(VoteDetails$Term > 1989 & VoteDetails$Term <2000, 1990, 2000)))))
pdf(file="Plots/TimeTrendsInTopics3a.pdf",18,12)
par(mfrow=c(4,6), mar=c(1,1,3,1), oma=c(4,4,1,1))
for (k in 1:Ntopics){
  beanplot(LDAIRT.out$topic.chain[k,,1] ~ VoteDecade
          ,main='',xlab="Year",ylab="Case Fraction"#,xaxs='i',yaxs='i'
          ,do.out=FALSE, do.conf=FALSE, outline=TRUE)
  title(main=list(LDATopicLabels[k], cex=1.6))
  #   polygon(c(1948,1953,1953,1948),c(0,0,.2,.2), col='grey',border=NA)
  #   polygon(c(1969,1986,1986,1969),c(0,0,.2,.2), col='grey',border=NA)
}
mtext(side=1,'Year',outer=TRUE, line=2)
mtext(side=2,'Case Faction',outer=TRUE, line=2)
dev.off()

pdf(file="Plots/TimeTrendsInTopics3.pdf",18,12)
par(mfrow=c(4,6), mar=c(1,1,3,1), oma=c(4,4,1,1))
for (k in 1:Ntopics){
  UniqueTerms <- unique(VoteDetails$Term)
  for(t in 1:length(unique(VoteDetails$Term))){
    tempKeepers <- which(VoteDetails$Term == UniqueTerms[t])
    tempProportions <- LDAIRT.out$topic.chain[k,tempKeepers,1]
    
    
    
  }
	plot(loess.out$x,loess.out$fitted,type="n",ylim=c(0,0.2),col="black",lwd=2,main='',xlab="Year",ylab="Mean Case Fraction",xaxs='i',yaxs='i')
	title(main=list(LDATopicLabels[k], cex=1.6))
	polygon(c(1948,1953,1953,1948),c(0,0,.2,.2), col='grey',border=NA)
	polygon(c(1969,1986,1986,1969),c(0,0,.2,.2), col='grey',border=NA)
	}
mtext(side=1,'Year',outer=TRUE, line=2)
mtext(side=2,'Case Faction',outer=TRUE, line=2)
dev.off()


pdf(file="Plots/TimeTrendsInTopics2S.pdf",9,3)
par(mfrow=c(1,3), mar=c(1,1,3,1), oma=c(4,4,1,1))
for (k in IssueSubsetForPresentation){
	loess.out <- loess(LDAIRT.out$topic.chain[k,,1]~VoteDetails$Term,span=0.5)
	plot(loess.out$x,loess.out$fitted,type="n",ylim=c(0,0.2),col="black",lwd=2,main='',xlab="Year",ylab="Mean Case Fraction",xaxs='i',yaxs='i')
	title(main=list(LDATopicLabels[k], cex=1.6))
	polygon(c(1948,1953,1953,1948),c(0,0,.2,.2), col='grey',border=NA)
	polygon(c(1969,1986,1986,1969),c(0,0,.2,.2), col='grey',border=NA)
	lines(loess.out$x,loess.out$fitted, lwd=2)
	}
mtext(side=1,'Year',outer=TRUE, line=2)
mtext(side=2,'Mean Case Faction',outer=TRUE, line=2)
dev.off()

for (k in 1:Ntopics){
	pdf(file=paste("Plots/TimeTrendsInTopics/",k,".pdf",sep=""),4,4)
	par(mfrow=c(4,6),mar=c(5.1,4.1,4.1,2.1))
	loess.out <- loess(LDAIRT.out$topic.chain[k,,1]~VoteDetails$Term,span=0.5)
	plot(loess.out$x,loess.out$fitted,type="l",ylim=c(0,0.2),col="black",lwd=2,main=LDATopicLabels[k],xlab="Year",ylab="Mean Case Fraction")
	dev.off()	
}

for (k in 1:Ntopics){
	pdf(file=paste("Plots/TimeTrendsInTopics/",k,"2.pdf",sep=""),4,4)
	par(mfrow=c(4,6),mar=c(3.1,4.1,3.1,2.1))
	loess.out <- loess(LDAIRT.out$topic.chain[k,,1]~VoteDetails$Term,span=0.5)
	plot(loess.out$x,loess.out$fitted,type="n",ylim=c(0,0.2),col="black",lwd=2,main=LDATopicLabels[k],xlab="Year",ylab="Mean Case Fraction", xaxs='i',yaxs='i')
	polygon(c(1948,1953,1953,1948),c(0,0,.2,.2), col='grey',border=NA)
	polygon(c(1969,1986,1986,1969),c(0,0,.2,.2), col='grey',border=NA)
	lines(loess.out$x,loess.out$fitted,col="black",lwd=2)
	box()
	dev.off()	
}
	
## Plot of Justice Positions in Each Issue ##	

pdf(file="Plots/JusticeLocationsInTopics.pdf",18,12)	
par(mfrow=c(4,6),cex.axis=0.75,mar=c(2.1,4.1,2.1,2.1))
sortindex <- order(colMeans(theta.est * rowSums(lambda.est)))
for (k in 1:Ntopics){
	plot(0,0,type="n",xlim=c(-4.5,4.5),ylim=c(0,Nvoters+1),xlab="",ylab="",main='',axes=FALSE)
	title(main=list(LDATopicLabels[k], cex=1.6))
	# axis(1)
	axis(2,at=1:Nvoters,labels=colnames(VoteMatrix)[sortindex],las=1)
	for (i in 1:Nvoters){
		lines(c(-10,10),c(i,i),col="lightgrey")
		points(theta.est[k,sortindex[i]],i,pch=16,cex=0.75)		
		lines(c(theta.est[k,sortindex[i]]-theta.se[k,sortindex[i]],theta.est[k,sortindex[i]]+theta.se[k,sortindex[i]]),c(i,i))

	
	}
}		
dev.off()	

pdf(file="Plots/JusticeLocationsInTopicsS.pdf",9,6)	
par(mfrow=c(1,3),cex.axis=1.25,mar=c(2.1,7.5,2.1,2.1))
sortindex <- order(colMeans(theta.est * rowSums(lambda.est)))
for (k in IssueSubsetForPresentation){
	plot(0,0,type="n",xlim=c(-4.5,4.5),ylim=c(0,Nvoters+1),xlab="",ylab="",main='',axes=FALSE)
	title(main=list(LDATopicLabels[k], cex=1.6))
	# axis(1)
	axis(2,at=1:Nvoters,labels=colnames(VoteMatrix)[sortindex],las=1)
	for (i in 1:Nvoters){
		lines(c(-10,10),c(i,i),col="lightgrey")
		points(theta.est[k,sortindex[i]],i,pch=16,cex=1)		
		lines(c(theta.est[k,sortindex[i]]-theta.se[k,sortindex[i]],theta.est[k,sortindex[i]]+theta.se[k,sortindex[i]]),c(i,i),lwd=2)

	
	}
}		
dev.off()	

for (k in 1:Ntopics){
	pdf(file=paste("Plots/JusticeLocationsInTopics/",k,".pdf",sep=""),5,5)	
	par(mfrow=c(1,1),cex.axis=0.7,mar=c(2.1,8.1,4.1,2.1))
	sortindex <- order(colMeans(theta.est * rowSums(lambda.est)))
	plot(0,0,type="n",xlim=c(-4.5,4.5),ylim=c(0,Nvoters+1),xlab="",ylab="",main=LDATopicLabels[k],axes=FALSE)
	# axis(1)
	axis(2,at=1:Nvoters,labels=colnames(VoteMatrix)[sortindex],las=1)
	for (i in 1:Nvoters){
		lines(c(-10,10),c(i,i),col="lightgrey")
		points(theta.est[k,sortindex[i]],i,pch=16,cex=0.75)		
		lines(c(theta.est[k,sortindex[i]]-theta.se[k,sortindex[i]],theta.est[k,sortindex[i]]+theta.se[k,sortindex[i]]),c(i,i))	
	}
	dev.off()	
}		

	
## Plot of Kennedy - O'Connor Posterior for Each Issue ##

pdf(file="Plots/KennedyOConnor.pdf",18,12)	
par(mfrow=c(4,6))
for (k in 1:Ntopics){
	plot(0,0,type="n",xlim=c(-2,2),ylim=c(0,2),xlab="Kennedy - O'Connor",ylab="Posterior Density",main=LDATopicLabels[k],axes=FALSE)	
	KennedyOConnor <- LDAIRT.out$ideal.chain[k,c(26),] -  LDAIRT.out$ideal.chain[k,c(24),]
	lines(density(KennedyOConnor))
	axis(1)	
	box()
	abline(v=0,col="grey")
	points(median(KennedyOConnor),0,pch=3)
	lines(quantile(KennedyOConnor,c(0.05,0.95)),c(0,0))
	text(1.5,1.75,format(mean(KennedyOConnor > 0),digits=3,nsmall=3))
}		
dev.off()

LDATopicLabelsMatrix <- matrix("",nrow=3,ncol=Ntopics)
for (k in 1:Ntopics) LDATopicLabelsMatrix[1:3,k] <- WordSparseMatrix$dimnames$Terms[TopWordsByTopic[,k]][1:3]

pdf(file="Plots/KennedyOConnor2.pdf",18,12)	
par(mfrow=c(4,6), mar=c(1,1,3,1), oma=c(4,4,1,1))
for (k in 1:Ntopics){
	plot(0,0,type="n",xlim=c(-2,2),ylim=c(0,2),xlab="",ylab="",main='',axes=FALSE)	
	title(main=list(LDATopicLabels[k], cex=1.6))
	KennedyOConnor <- LDAIRT.out$ideal.chain[k,c(26),] -  LDAIRT.out$ideal.chain[k,c(24),]
	lines(density(KennedyOConnor))
	axis(1)	
	box()
	abline(v=0,col="grey")
	points(median(KennedyOConnor),0,pch=3)
	lines(quantile(KennedyOConnor,c(0.05,0.95)),c(0,0))
	text(1.5,1.75,format(mean(KennedyOConnor > 0),digits=3,nsmall=3))
}
	mtext(side=1,"Kennedy-'O'Connor", outer=TRUE, line=2)
	mtext(side=2,"Posterior Density", outer=TRUE, line=2)		
dev.off()

KennedyOConnor <- matrix(NA,nrow=Ntopics,ncol=dim(LDAIRT.out$ideal.chain)[3])
KennedyOConnorMedian <- NULL
for (k in 1:Ntopics){
	KennedyOConnor[k,] <- LDAIRT.out$ideal.chain[k,c(26),] -  LDAIRT.out$ideal.chain[k,c(24),]
	KennedyOConnorMedian[k] <- median(KennedyOConnor[k,])
}
KennedyOConnorOrder <- order(KennedyOConnorMedian, decreasing=FALSE)
pdf('Plots/KennedyOConnorSingleFigure.pdf',6,8)
par(mar=c(5,13,3,1))
plot(0,0,type="n",xlim=c(-3,3),ylim=c(0,Ntopics),xlab="",ylab="",main='',axes=FALSE)
# axis(1)
axis(2,at=1:Ntopics,labels=LDATopicLabels[KennedyOConnorOrder],las=1)
for (i in 1:Ntopics){
	lines(c(-10,10),c(i,i),col="lightgrey")
	points(KennedyOConnorMedian[KennedyOConnorOrder[i]],i,pch=16,cex=1)		
	lines(c(quantile(KennedyOConnor[KennedyOConnorOrder[i],], probs=c(.05)),quantile(KennedyOConnor[KennedyOConnorOrder[i],], probs=c(.95))), c(i,i),lwd=2)
}
lines(c(0,0), c(1,Ntopics), lty=5, col='grey')
axis(side=1)
mtext(side=3, line = 0, at=-2, 'Kennedy More Liberal', cex=.7)
mtext(side=3, line = 0, at=2, 'Kennedy More Conservative', cex=.7)
mtext(side=1,line=2,'Posterior Difference Between Kennedy and O\'Connor', cex=.7)
dev.off()

for (k in 1:Ntopics){
	pdf(file=paste("Plots/KennedyOConnor/",k,".pdf",sep=""),4,4)	
	par(mfrow=c(1,1))
	plot(0,0,type="n",xlim=c(-2,2),ylim=c(0,2),xlab="Kennedy - O'Connor",ylab="Posterior Density",main=LDATopicLabels[k],axes=FALSE)	
	KennedyOConnor <- LDAIRT.out$ideal.chain[k,c(26),] -  LDAIRT.out$ideal.chain[k,c(24),]
	lines(density(KennedyOConnor))
	axis(1)	
	box()
	abline(v=0,col="grey")
	points(median(KennedyOConnor),0,pch=3)
	lines(quantile(KennedyOConnor,c(0.05,0.95)),c(0,0))
	text(1.5,1.75,format(mean(KennedyOConnor > 0),digits=3,nsmall=3))
	dev.off()
}		


### make a plot of median by topic by term
medianMatrix <- matrix(NA,nrow=24, ncol=LDAIRT.summary$Nvotes)
for(i in 1:24){
	for(j in 1:LDAIRT.summary$Nvotes){
		medianMatrix[i,j] <- median(theta.est[i,which(!is.na(VoteMatrix[j,]))])
	}
}

courtComposition <- read.csv('compositionByTerm.csv', header=TRUE, stringsAsFactors=FALSE, na.strings='NA')

termVector <- unique(VoteDetails$Term)
medianMatrix <- matrix(NA, nrow=24, ncol=length(termVector))
for(i in 1:24){
	for(j in 1:length(termVector)){
		medianMatrix[i,j] <- median(LDAIRT.summary$theta.est[i,which(colnames(LDAIRT.summary$theta.est) %in% courtComposition[j,2:10])
		])
	}
}

pdf(file="Plots/MedianByTopic.pdf",18,12)	
par(mfrow=c(4,6), mar=c(1,1,3,1), oma=c(4,4,1,1))
for (k in 1:Ntopics){
	plot(termVector,medianMatrix[k,],type="n",ylim=c(-2,2),col="black",lwd=2,main='',xlab="",ylab="",xaxs='i',yaxs='i')
	title(main=list(LDATopicLabels[k], cex=1.6))
	polygon(c(1948,1953,1953,1948),c(-2,-2,2,2), col='grey',border=NA)
	polygon(c(1969,1986,1986,1969),c(-2,-2,2,2), col='grey',border=NA)
	lines(termVector,medianMatrix[k,], lwd=2)
	}
mtext(side=1,'Year',outer=TRUE, line=2)
mtext(side=2,'Median Justice Conservatism',outer=TRUE, line=2)
dev.off()

pdf(file="Plots/MedianByTopicS.pdf",9,3)	
par(mfrow=c(1,3))
for (k in IssueSubsetForPresentation){
	plot(termVector,medianMatrix[k,],type="n",ylim=c(-2,2),col="black",lwd=2,main='',xlab="",ylab="",xaxs='i',yaxs='i')
	title(main=list(LDATopicLabels[k], cex=1.6))
	polygon(c(1948,1953,1953,1948),c(-2,-2,2,2), col='grey',border=NA)
	polygon(c(1969,1986,1986,1969),c(-2,-2,2,2), col='grey',border=NA)
	lines(termVector,medianMatrix[k,], lwd=2)
	}
mtext(side=1,'Year',outer=TRUE, line=2)
mtext(side=2,'Median Justice Conservatism',outer=TRUE, line=-1.5)
dev.off()




# # Kyllo still doesn't work well!
# pdf('Plots/1DIdealPointPlotC.pdf',8,4)
# par(mfrow=c(1,1),mar=c(4.1,4.1,2.1,2.1))
# SpecificCaseVotes <- VoteMatrix[4141,22:30]
# SpecificCaseDiscrim <- LDAIRT.out$discrim.chain[4141,,]
# CaseSpecificIdeal <- lambda.est[,4141] %*% theta.est[,22:30] 
# plot(CaseSpecificIdeal,SpecificCaseVotes,type="n",axes=FALSE,xlab="1D Justice Ideal Point",ylab="Justice's Dispositional Vote",main="U.S. Supreme Court Justice Votes\n in Kyllo v. US (2001) by 1D Ideal Point",ylim=c(-0.1,1.5),xlim=c(-0.5,3.5))
# xseq <- seq(-0.5,3.5,0.05)
# pseq <- 1-pnorm(xseq*SpecificCaseDiscrim[2] + SpecificCaseDiscrim[1]) 
# lines(xseq,pseq,col="grey")
# points(CaseSpecificIdeal,SpecificCaseVotes,pch=16)
# text(CaseSpecificIdeal,SpecificCaseVotes-0.15*(SpecificCaseVotes-0.5),colnames(CaseSpecificIdeal)[22:30],pos=4-2*SpecificCaseVotes,srt=60,offset=-0.2,cex=0.8)
# axis(1,labels=FALSE)
# axis(2,at=c(0,1),labels=c("US","Kyllo"))
# dev.off()




### 2D Model ###


## Load Best MCMC Chain and Create Summaries ## 

load(file="OptimalTopicCountSearch/SavedLDAIRT-2-5000-5000.RData")
LDAIRT.summary <- summary.LDAIRT(LDAIRT.out)
attach(LDAIRT.summary)
theta.se <- apply(LDAIRT.out$ideal.chain,c(1,2),sd)

SpaethCategoryNames <- c("Criminal Procedure","Civil Rights","First Amendment","Due Process","Privacy","Attorneys","Unions","Economic Activity","Judicial Power","Federalism","Interstate Relations","Federal Taxation","Miscellaneous")
LDATopicLabels <- rep("",Ntopics)
for (k in 1:Ntopics) LDATopicLabels[k] <- paste((WordSparseMatrix$dimnames$Terms[TopWordsByTopic[,k]])[1:5],collapse=", ")
LDATopicTable <- matrix("",10,Ntopics)
for (k in 1:Ntopics) LDATopicTable[,k] <- WordSparseMatrix$dimnames$Terms[TopWordsByTopic[1:10,k]]

pdf('Plots/2DIdealPointPlot.pdf',16,8)
par(mfrow=c(1,2))
plot(theta.est[1,],theta.est[2,],xlim=c(-2,4.1),ylim=c(-2,4.1), type='n', xlab=paste('Dimension 1:',LDATopicLabels[1]), ylab=paste('Dimension 2:',LDATopicLabels[2])
	,main='All Justices (except Douglas)')
text(theta.est[1,], theta.est[2,], labels=colnames(theta.est), cex=.5)
for(i in 1:dim(LDAIRT.out$ideal.chain)[2]){
	dataEllipse(LDAIRT.out$ideal.chain[1,i,],LDAIRT.out$ideal.chain[2,i,], levels=c(0.9),plot.points=FALSE,col="light grey",lwd=0.5,center.cex=1,center.pch=0)
}
plot(theta.est[1,],theta.est[2,],xlim=c(-1.1,1.1),ylim=c(-2,2), type='n', xlab=paste('Dimension 1:',LDATopicLabels[1]), ylab=paste('Dimension 2:',LDATopicLabels[2])
	,main='Moderate Justices')
text(theta.est[1,], theta.est[2,], labels=colnames(theta.est), cex=.5)
for(i in 1:dim(LDAIRT.out$ideal.chain)[2]){
	dataEllipse(LDAIRT.out$ideal.chain[1,i,],LDAIRT.out$ideal.chain[2,i,], levels=c(0.9),plot.points=FALSE,col="light grey",lwd=0.5,center.cex=1,center.pch=0)
}
dev.off()




### 1D Model ###


## Load Best Chain and Create Summaries ## 

load(file="OptimalTopicCountSearch/SavedLDAIRT-1-2000-2000.RData")
theta.est <- apply(LDAIRT.out$ideal.chain,c(1,2),mean)
theta.se <- apply(LDAIRT.out$ideal.chain,c(1,2),sd)
theta.q <- apply(LDAIRT.out$ideal.chain,c(1,2),quantile,c(0.025,0.975))


pdf('Plots/1DIdealPointPlotA.pdf',8,4)
par(mfrow=c(1,1),mar=c(4.1,4.1,2.1,2.1))
SpecificCaseVotes <- VoteMatrix[4123,22:30]
SpecificCaseDiscrim <- LDAIRT.out$discrim.chain[4123,,]
plot(-theta.est[22:30],SpecificCaseVotes,type="n",axes=FALSE,xlab="1D Justice Ideal Point",ylab="Justice's Dispositional Vote",main="U.S. Supreme Court Justice Votes\n in Bush v. Gore (2000) by 1D Ideal Point",ylim=c(-0.1,1.5),xlim=c(-0.5,3.5))
xseq <- seq(-0.5,3.5,0.05)
pseq <- 1-pnorm(xseq*SpecificCaseDiscrim[2] + SpecificCaseDiscrim[1]) 
lines(xseq,pseq,col="grey")
points(-theta.est[22:30],SpecificCaseVotes,pch=16)
text(-theta.est[22:30],SpecificCaseVotes-0.15*(SpecificCaseVotes-0.5),colnames(theta.est)[22:30],pos=4-2*SpecificCaseVotes,srt=60,offset=-0.2,cex=0.8)
axis(1,labels=FALSE)
axis(2,at=c(0,1),labels=c("Gore","Bush"))
dev.off()


pdf('Plots/1DIdealPointPlotB.pdf',8,4)
par(mfrow=c(1,1),mar=c(4.1,4.1,2.1,2.1))
SpecificCaseVotes <- VoteMatrix[4141,22:30]
SpecificCaseDiscrim <- LDAIRT.out$discrim.chain[4141,,]
plot(-theta.est[22:30],SpecificCaseVotes,type="n",axes=FALSE,xlab="1D Justice Ideal Point",ylab="Justice's Dispositional Vote",main="U.S. Supreme Court Justice Votes\n in Kyllo v. US (2001) by 1D Ideal Point",ylim=c(-0.1,1.5),xlim=c(-0.5,3.5))
xseq <- seq(-0.5,3.5,0.05)
pseq <- 1-pnorm(xseq*SpecificCaseDiscrim[2] + SpecificCaseDiscrim[1]) 
lines(xseq,pseq,col="grey")
points(-theta.est[22:30],SpecificCaseVotes,pch=16)
text(-theta.est[22:30],SpecificCaseVotes-0.15*(SpecificCaseVotes-0.5),colnames(theta.est)[22:30],pos=4-2*SpecificCaseVotes,srt=60,offset=-0.2,cex=0.8)
axis(1,labels=FALSE)
axis(2,at=c(0,1),labels=c("US","Kyllo"))
dev.off()




pdf('Plots/1DIdealPointPlotC.pdf',8,4)
par(mfrow=c(1,1),mar=c(4.1,4.1,2.1,2.1))

whichCase <- 1938
whichCaseName <- "Miller v. California (1973)"
whichJustices <- is.na(VoteMatrix[whichCase,]) == FALSE

load(file="OptimalTopicCountSearch/SavedLDAIRT-1-2000-2000.RData")
theta.est <- apply(LDAIRT.out$ideal.chain,c(1,2),mean)
theta.se <- apply(LDAIRT.out$ideal.chain,c(1,2),sd)
theta.q <- apply(LDAIRT.out$ideal.chain,c(1,2),quantile,c(0.025,0.975))

DiscrimDiff <- -LDAIRT.out$discrim.chain

SpecificCaseVotes <- VoteMatrix[whichCase,]
SpecificCaseDiscrim <- LDAIRT.out$discrim.chain[whichCase,,]
plot(-theta.est,SpecificCaseVotes,type="n",axes=FALSE,xlab="1D Justice Ideal Point",ylab="Justice's Dispositional Vote",main="U.S. Supreme Court Justice Votes\n in Miller v. California (1973) by 1D Ideal Point",ylim=c(-0.1,1.5),xlim=c(-8,3.5))
xseq <- seq(-7.5,3.5,0.05)
pseq <- 1-pnorm(xseq*SpecificCaseDiscrim[2] + SpecificCaseDiscrim[1]) 
lines(xseq,pseq,col="grey")
points(-theta.est,SpecificCaseVotes,pch=16)
text(-theta.est[whichJustices],SpecificCaseVotes[whichJustices]-0.15*(SpecificCaseVotes[whichJustices]-0.5),colnames(theta.est)[whichJustices],pos=4-2*SpecificCaseVotes[whichJustices],srt=60,offset=-0.2,cex=0.8)
axis(1,labels=FALSE)
axis(2,at=c(0,1),labels=c("Miller","California"))

dev.off()


pdf('Plots/1DIdealPointPlotC2.pdf',8,8)

par(mfrow=c(2,1),mar=c(4.1,4.1,2.1,2.1))


whichCase <- 1938
whichCaseName <- "Miller v. California (1973)"
whichJustices <- is.na(VoteMatrix[whichCase,]) == FALSE

load(file="OptimalTopicCountSearch/SavedLDAIRT-1-2000-2000.RData")
theta.est <- apply(LDAIRT.out$ideal.chain,c(1,2),mean)
theta.se <- apply(LDAIRT.out$ideal.chain,c(1,2),sd)
theta.q <- apply(LDAIRT.out$ideal.chain,c(1,2),quantile,c(0.025,0.975))

DiscrimDiff <- -LDAIRT.out$discrim.chain

SpecificCaseVotes <- VoteMatrix[whichCase,]
SpecificCaseDiscrim <- LDAIRT.out$discrim.chain[whichCase,,]
plot(-theta.est,SpecificCaseVotes,type="n",axes=FALSE,xlab="1D Justice Ideal Point",ylab="Justice's Dispositional Vote",main=whichCaseName,ylim=c(-0.1,1.5),xlim=c(-8,3.5))
xseq <- seq(-7.5,3.5,0.05)
pseq <- 1-pnorm(xseq*SpecificCaseDiscrim[2] + SpecificCaseDiscrim[1]) 
lines(xseq,pseq,col="grey")
points(-theta.est,SpecificCaseVotes,pch=16)
text(-theta.est[whichJustices],SpecificCaseVotes[whichJustices]-0.15*(SpecificCaseVotes[whichJustices]-0.5),colnames(theta.est)[whichJustices],pos=4-2*SpecificCaseVotes[whichJustices],srt=60,offset=-0.2,cex=0.8)
axis(1,labels=FALSE)
axis(2,at=c(0,1),labels=c("Miller","California"))

load(file="OptimalTopicCountSearch/SavedLDAIRT-24-5000-5000.RData")
LDAIRT.summary <- summary.LDAIRT(LDAIRT.out)
theta.est <- LDAIRT.out$topic.chain[,whichCase,1] %*% apply(LDAIRT.out$ideal.chain,c(1,2),mean)

DiscrimDiff <- LDAIRT.out$discrim.chain - DiscrimDiff

SpecificCaseVotes <- VoteMatrix[whichCase,]
SpecificCaseDiscrim <- LDAIRT.out$discrim.chain[whichCase,,]
plot(theta.est,SpecificCaseVotes,type="n",axes=FALSE,xlab="Case-Specific 1D Justice Ideal Point",ylab="Justice's Dispositional Vote",main="",ylim=c(-0.1,1.5),xlim=c(-8,3.5))
xseq <- seq(-7.5,3.5,0.05)
pseq <- 1-pnorm(-xseq*SpecificCaseDiscrim[2] + SpecificCaseDiscrim[1]) 
lines(xseq,pseq,col="grey")
points(theta.est,SpecificCaseVotes,pch=16)
text(theta.est[whichJustices]+c(0,0,0,0,0,0,-0.2,0,0),SpecificCaseVotes[whichJustices]-0.15*(SpecificCaseVotes[whichJustices]-0.5),colnames(theta.est)[whichJustices],pos=4-2*SpecificCaseVotes[whichJustices],srt=60,offset=-0.2,cex=0.8)
axis(1,labels=FALSE)
axis(2,at=c(0,1),labels=c("Miller","California"))


dev.off()


### analysis of vote dispositions before and after o'connor's departure
load('SpaethData/SCDB_2011_03_caseCentered_Citation.Rdata')
spaethData <- SCDB_2011_03_caseCentered_Citation
rm(SCDB_2011_03_caseCentered_Citation)	

term <- NULL
for(i in 1:length(colnames(LDAIRT.summary$lambda.est))){
	term[i] <- spaethData$term[which(spaethData$usCite == colnames(LDAIRT.summary$lambda.est)[i])]
}
	
	
	